# include "scene.h"
# include <math.h>
# include <stdlib.h>

/* This function returns inlet velocity profile and its time derivatives. 
   w(r,t) = W_MAX*[1 - (r/r_NOZZLE)^2]*[1 - e^(-b.t)]. */

extern double inlet_velocity(double r, double t, short index)
     {
      double b, q;
/*
       double w;

       r = r*1000; //convert r into mm.
       w = (-0.1762*r*r*r*r + 1.3982*r*r*r - 6.6613*r*r - 0.8037*r + 102.16)/1000; //Vel in "m/sec".
       return w;
*/
    /* For negative time, return only radial function of the velocity. */
//        if(t < -1*NEG_DELTA) return W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE));

       b = -log(1-wFRAC)/tFRAC;

    /* Return integration of the velocity from time "0 ~ t". */
       if(index == -1) q = W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*(t + (exp(-b*t) - 1)/b);

       if(index > 3)
          {
           printf("ERROR:: In 'inlet_velocity()': 'index = %d'. However, possible values of 'index': 0, 1, 2, or 3.\
                  'index = n' will return n-th time-derivative of the function.\n", index);
           exit(-1);
          }

    /* return w(r,t). */
       if(index == 0) q = W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*(1 - exp(-b*t));

    /* return dw(r,t)/dt. */
       if(index == 1) q = b*W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*exp(-b*t);

    /* return d^2w(r,t)/dt^2. */
       if(index == 2) q = -b*b*W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*exp(-b*t);

    /* return d^3w(r,t)/dt^3. */
       if(index == 3) q = b*b*b*W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*exp(-b*t);

    /* return d^3w(r,t)/dt^4. */
//        if(index == 4) q = -b*b*b*b*W_MAX*(1 - (r/r_NOZZLE)*(r/r_NOZZLE))*exp(-b*t);

       return q;
     }
